function [A,B,CF]=csr(A,B,CF,nf,nnf,nleads);

% CSR.M: classical state reduction program
% inputs: A,B,CF,nf,nleads
% outputs: A,B,CF

% The following program solves a classical state reduduction
% problem for the dynamic system of the form:
%   A Ey(t+1|t) = B y(t) + C(F) Ex(t|t)    
% with the posited structure

%  |0  0   0  | |Eo(t+1|t)|   |I   Bof Bod| |o(t)|    |Co(F)|
%  |          | |         |   |           | |    |  + |     |  
%  |0  0   0  | |Ef(t+1|t)| = |0   Bff Bfd| |f(t)|    |Cf(F)| Ex(t|t) 
%  |          | |         |   |           | |    |  + |     |     
%  |0  Adf Add| |Ed(t+1|t)|   |0   Bdf Bdd| |d(t)|    |Cd(F)| 

% In these expressions, o(t) are previously isolated flows ("old flows")
% f(t) are new flows d(t) are dynamic variables.  It is assumed that
% Bff is a nonsingular (nf by nf) matrix; o(t) is the first nf-nnf
% elements of y(t), f(t) is the next nnf elements and 
% d(t) is the last nd=ny-nf elements of y(t). 

ny=rows(A);
nd=ny-nf;
nof=nf-nnf;

% We compute the state reduction using our knowledge of this special
% structure rather than the equivalent transformation T(F) = G*F+H 
% described in the manuscript, so as to eliminate large matrix
% multiplications.  However, the net result is that the new system
% is in the form:

%  | 0  0 | |Ef(t+1|t)|   |I     NU | |f(t)|    |CCf(F)|
%  |      | |         | = |         | |    |  + |      |  Ex(t|t)    
%  |0  AA | |Ed(t+1|t)|   |0     BB | |d(t)|    |CCd(F)| 

% with f(t) now containing all the o(t) and f(t) elements above.
% That is, the transformed system has specified flow variables and 
% state variables.  We do not require that AA is invertible, so that the
% current program can be used in a singular systems context.

% Define matrices (locally) for programming convenience

Bff=B(nof+1:nf,nof+1:nf);
Bfd=B(nof+1:nf,nf+1:ny);
CFf=CF(nof+1:nf,:);

% Step 1: normalize the new flow equations.

B(nof+1:nf,nof+1:nf)=eye(nnf);
B(nof+1:nf,nf+1:ny)=Bff\Bfd;
CF(nof+1:nf,:)=Bff\CFf;

Bfdnew=B(nof+1:nf,nf+1:ny);
CFfnew=CF(nof+1:nf,:);

% This step puts the system in the form

%  |0  0   0  | |Eo(t+1|t)|   |I  Bof Bod| |o(t)|    |Co(F)|
%  |          | |         |   |          | |    |  + |     |  
%  |0  0   0  | |Ef(t+1|t)| = |0  I   Bfd| |f(t)|    |Cf(F)| Ex(t|t) 
%  |          | |         |   |          | |    |  + |     |     
%  |0  Adf Add| |Ed(t+1|t)|   |0  Bdf Bdd| |d(t)|    |Cd(F)| 

% with suitable redefinitions of matrices. We call these "new"

% Step 2: eliminate the dependence of o(t) on f(t), 
%         subtract Bof times second block to first block


Bod=B(1:nof,nf+1:ny);
Bof=B(1:nof,nof+1:nf);

B(1:nof,nof+1:nf)=zeros(nof,nnf);
B(1:nof,nf+1:ny)=Bod-Bof*Bfdnew;

% CFo(new)=  CFo(old) - Bof*CFf(new)

CF(1:nof,:)=CF(1:nof,:)-Bof*CFfnew;

% This puts the  system in the form:
%  |0  0   0  | |Eo(t+1|t)|   |I  0   Bod| |o(t)|    |Co(F)|
%  |          | |         |   |          | |    |  + |     |  
%  |0  0   0  | |Ef(t+1|t)| = |0  I   Bfd| |f(t)|    |Cf(F)| Ex(t|t) 
%  |          | |         |   |          | |    |  + |     |     
%  |0  Adf Add| |Ed(t+1|t)|   |0  Bdf Bdd| |d(t)|    |Cd(F)| 

% with a new definition of Bod.

% Step 3: eliminate the dependence of d(t+1) on f(t+1)
% and f(t).

% Add Adf times the "updated" second block (as previously transformed)
% to the third block

Add=A(nf+1:ny,nf+1:ny);
Adf=A(nf+1:ny,nof+1:nf);

A(nf+1:ny,nf+1:ny)=Add-Adf*Bfdnew;
A(nf+1:ny,nof+1:nf)=zeros(nd,nnf);

% Subtract Bdf times the second block (as previously transformed)
% to the third block

Bdd=B(nf+1:ny,nf+1:ny);
Bdf=B(nf+1:ny,nof+1:nf);

B(nf+1:ny,nf+1:ny)=Bdd-Bdf*Bfdnew;
B(nf+1:ny,nof+1:nf)=zeros(nd,nnf);

% make the new C(F) polynomial reflect these two transformations
% which are in the form of a lead polynomial

nx=cols(CF)/(nleads+1);

CFd=[CF(nf+1:ny,:) zeros(nd,nx)]+...
              mconv([-Bdf Adf],CFfnew,1,nleads);

CF=[CF(1:nf,:) zeros(nf,nx)
    CFd];


